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We propose a modified Gaussian ansatz to study binary condensates, trapped in harmonic and 
optical lattice potentials, both in miscible and immiscible domains. The ansatz is an apt one as it 
leads to the smooth transition from miscible to immiscible domains without any a priori assumptions. 
In optical lattice potentials, we analyze the squeezing of the density profiles due to the increase in 
the depth of the optical lattice potential. For this we develop a model with three potential wells, 
and define the relationship between the lattice depth and profile of the condensate. 
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I. INTRODUCTION 

After the successful experimental realization of two- 
species Bose-Einstein condensate (TBEC), consisting of 
different hyperfine spin states of 87 Rb P~H3], experimental 
as well as theoretical investigations in this field have come 
a long way. TBECs of different atomic species ( 41 K and 
87 Rb) 4j and of different isotopes of same atomic species 
( 85 Rb and 87 Rb) [5] have been experimentally observed. 
Phase separation, a typical feature of two component 
Bose-Einstein condensates, has been observed unambigu- 
ously [5]. A lot of static as well as dynamical properties 
of the TBECs have been analyzed in great detail in recent 
years. These include the ground state geometry [Btil4|. 
modulational instability [T5lll7| . Rayleigh- Taylor insta- 
bility Ql], [IS], Kelvin-Helmholtz instability gfr], etc.. 

The ground state geometry of the TBEC has been 
studied, semi-analytically, using Thomas-Fermi (TF) ap- 
proximation [BH5] . These works did not take into account 
the contribution of the interface energy explicitly. The 
interface energy was taken into account in some other 
works [TUHH]; however, the interface energy correction 
incorporated in these works was not good enough to re- 
produce the experimentally observed ground state struc- 
tures with planar interface [5] . A more accurate analytic 
approximation to account for interface energy was sug- 
gested by Schaeybroeck et al. [13] . and a recent work 
|14j conclusively proved that using this analytic approx- 
imation for interface energy in TF regime, planar and 
cylindrical geometries emerge as the ground state struc- 
tures in cigar and pan-cake shaped trap potentials. The 
common salient feature of all these works is the use of TF 
approximation, which is a good approximation for large 
number of atoms of each species. In fact if A, a, and 
a osc are respectively, number of atoms of the component 
species, s-wave scattering length, and oscillator length, 
TF approximation is valid provided Na/a osc 3> 1. Ob- 
viously, this condition is not satisfied for attractive con- 
densates. For small number of atoms as well, say of the 
order of a few hundreds, and for very weakly interacting 
condensates (a/a osc <C 1), the contribution of the kinetic 
energy to the total energy is significant and can not be 



merely treated as a correction to total energy as is done 
in TF based approaches. For such TBECs, TF approxi- 
mation is not a good approximation, and hence can not 
be relied upon to determine the ground state structure 
of the TBEC. In present work, we analyze the ground 
sate properties of the TBECs using suitable ansatz in 
both miscible and immiscible regimes. The Gaussian na- 
ture of the ansatz makes our approach better equipped to 
study the stationary state properties of the of very weakly 
interacting (Na/a osc ~ 1) TBECs and also those with 
attractive interactions. With the advent of Feshbach res- 
onances [21] , it is experimentally possible to tune scatter- 
ing lengths to reach the very weakly interacting regime or 
even non-interacting regime [2"2T - [2"4"] . Magnetic Feshbach 
resonances can tune only one scattering length indepen- 
dently, whereas optical Feshbach resonances [25] open 
up the possibility of tuning different scattering length 
in a multicomponent system independently. With the 
experimental observation of optical Feshbach resonances 
in 172 Yb[53] and 174 Yb [57], experimental realization of 
weakly interacting regime in binary condensates appears 
a distinct possibility. The present approach, thus, sup- 
plements TF based semi-analytic schemes to determine 
the ground state geometries of the TBECs. 



Keeping a pace with the studies on TBECs, has been 
the realization of various condensed matter phenom- 
ena like ac Josephson effect, Bloch oscillations, Landau- 
Zcner tunneling, etc. in single species Bose-Einstein 
condensates (BECs) trapped in optical lattices [28ll31j . 
Superfluid-Mott insulator transition has been also ob- 
served with BECs in optical lattice potentials [35]. In 
mean field approximation, we also study the ground sate 
geometry of TBECs trapped in optical lattice potentials 
using discrete nonlinear Schrodinger equation (DNLSE) 
[33] , We also compare our semi-analytic results with the 
numerical solution of coupled Gross-Pitaevskii equations 
and find a very good agreement especially in very weakly 
interacting regime. 
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II. TBECS IN AXISYMMETRIC TRAPS 

In this section, we provide a general variational scheme 
to study the stationary state geometry of TBECs in ax- 
isymmetric trap potentials 



Vi(r,z) 



rriiUJ 



(r 2 + afz 2 ) 



(1) 



where i = 1,2 is the species index, uj is the radial trap 
frequency for two components, and Ui are the anisotropy 
parameters. For simplicity of analysis, we consider trap 
potentials for the component species are identical a± = 
ai2 = a and mi = m 2 = m. The ground state of the 
TBEC is described by a set of coupled GP equations 



2m 



S7 2 + V i (r,z) + J2u ij \^ j 



(2) 



in mean held approximation, where i = 1, 2 is the species 
index. Here Uu = 4irh 2 au/mi, where m^ is the mass 
and an is the s-wave scattering length, is the intra- 
species interaction, [Zy = 2irh 2 a i j/m i j, where my — 
raiTtij I (mi + irij) is the reduced mass and ay is the inter- 
species scattering length, is the inter-species interaction, 
and fii is the chemical potential of the ith species. The 
energy of the TBEC is 



E = 



E 



2m 



V*i| 2 +Vi(r,z)* 4 2 



i7l2|^l| 2 |*2| 2 



dr. 



(3) 



To rewrite the energy in suitable units, define the oscil- 
lator length of the trapping potential 



(4) 



and consider tno as the unit of energy. We then divide 
the Eq.([3]) by Hlo and apply the transformations 



t = tlx), and E = 



E 



The transformed order parameter 



— z), 



(5) 



(6) 



and energy of the TBEC in scaled units is 



E 




-Vi(r,z)\<Pi\ 2 - 
NiN 2 Uia\<h.\ 2 \<h\ 2 \ , 



(7) 



where Uu = AiTaii/a osc and U\2 = ^a\2 /o, osc in the 
scaled units. For simplicity of notations, from here on 
we represent the transformed quantities without tilde. 
We obtain the coupled 2D GP equations, 



1^ 

2cV2 



i__a 

2r dr 



Vi(r,z) + y] Gjj \<pj\ 2 <t>i = iM<t>i, 

(8) 

when the energy functional £ = E — /XjiVj is varia- 
tionally minimized with </>* as parameters of variation, 



here Gi 



NiUu and G, 



N U l3 . 



The equations 



are then solved numerically. Another approach ideal for 
semi-analytic treatment is to adopt a predefined form of 
<t> with few variational parameters and minimize E. This 
is outlined for the 2D and ID in the next subsections. 



A. Variational ansatz 

As mentioned earlier, the ground state of TBEC can be 
either be miscible or immiscible depending on the inter- 
action parameters. An ansatz which describes the ground 
state of the TBEC well, both in the miscible and immis- 
cible domain, when an, ay > and 022 > an is 

9 2 (r,z) = [f + c(r 2 + a 2 z 2 )]e-^ +a2 ^^ 2d2 \ (9) 

with a, b, /, c, and d as the variational parameters. In 
the phase separated or immiscible domain, this ansatz 
is apt for the ellipsoidal interface geometry where the 
density distribution follows the equipotential surfaces of 
the trapping potentials. It is applicable, in particular, 
to weakly interacting TBECs, and this is precisely the 
underlying assumption for choosing the present ansatz. 
This ansatz is not suitable for analysis of planar and 
cylindrical geometries where density distributions do not 
follow equipotential surfaces. A similar ansatz, perhaps 
more general, is used in ref. [34] to examine symbiotic 
gap and semi gap solitons in TBECs trapped in opti- 
cal lattice potentials. The parameter c accounts for the 
flattening of the density profile of the second species as 
the intra-species non-linearity is increased. Moreover, for 
symmetric ground state geometries, c is a parameter re- 
lated to the overlap of two component wave functions. 
For an ideal case when all the non-linearities are small 
and equal, c is ~ 0, and hence the two species of the 
TBEC completely overlap. 

If N\ and N 2 are the number of particles of the two 
species, then from Eq.(|6| in scaled units 



dv\ ( j> i {r,z)\ 2 = l. 



(10) 



Since the number of atoms of each species is fixed, the 
normalization conditions are equivalent to two constraint 
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equations and reduce the number of variational parame- 
ters by two. After evaluating the integrals, the two con- 
straint equations are 

(-l) 2 /3 a i/3 



2^/3 / y / gg7 2 a - 2d 10 / 2 7r 3 r- 5 \ 



.(11) 



From Eqs.([7|T0|, the energy of the first species is 
a 2 bNnr^ 2 



so 



\/2a 2 & 2 Gii+4 + 6fe 4 + 2a 2 ) , (12) 



similarly, for the second species 

E 2 = dN2nV2 {l28 [4f 2 (2 + 3,/ 
2048a 1 



4 ' a 2 ) 



4cd 2 /(2 



+15d 4 + a 2 ) + c 2 d 4 (22 + 105d 4 + 11a 2 )] 
+y2d 2 G 22 (945c 4 d 8 + 1680c 3 d 6 / + 1440c 2 d 4 / 2 
+768cd 2 / 3 + 256/ 4 )} , (13) 

and the energy from the inter species interaction is 

2N 1 G 12 b 2 d 3 iT^ 2 



E\2 — 



[I56 4 c 2 d 4 + I2b 2 cd 2 f(b 2 



2048a(6 2 + d 2 ) 7 / 2 
+d 2 ) + 4f{b 2 + d 2 ) 2 ] . 

We can then define the energy per boson as 

_ Ei + E 2 + E\2 
6 ~ Wi +N 2 



(14) 



(15) 



This can now be minimized numerically to determine 
variational parameters and hence the ground state wave 
functions. The results of minimization with the param- 
eters satisfying Nia^ / a osc ~ 1 are shown in FigJT] along 
with the corresponding numerical results. 



It should be noted that the scenario of all the non- 
linearities to be equal (iVi = N 2 and an = a,j) is not 
equivalent to single component condensate with an — 
a%j — a$ as the s-wave scattering length and total num- 
ber of atoms equal to N\ + N 2 . This is due to fact that 
we are still treating the two components as two different 
species having order parameters 4>i( r i z ) an d 4>2( r , z )- ft 
means that experimentally, even if for two different hy- 
perfine states of an isotope above condition is satisfied, 
the system will be still a binary system having a pair of 
ground state wave functions 4>\(r, z) and 4> 2 {r, z) instead 
of single ground state wave function for single compo- 
nent BEC. Quantum mechanically, it means that for the 
system to behave as a single species BEC, the two com- 
ponents needs to be indistinguishable with the same wave 
function <p{r, z). 
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FIG. 1: The semi-analytic and numerical profiles for a binary 
condensate with iVi = N 2 = 10, 000, a = 0.8, a osc = 9.566 x 
10~ r m, an = 0.5ao, (122 = l.Oao, and 012 = 0.8ao. Starting 
from left, the upper panel shows the semi-analytic profiles for 
the first and the second component, exactly below them in 
lower panel are wave function profiles obtained by numerically 
solving Eq.|8]|. 



B. Quasi-ID condensates 

When the radial trapping frequency is much larger 
than axial trapping frequency (a <C 1), and the TBEC 
is in the weakly interacting regime a i N i \tp(z)\ 2 <C 1, the 
order parameter can be factorized into radial and axial 
parts 



<Ai(r, z) = £i(r)ipi(z), 



(16) 



where £i(r) is the normalized ground state of radial trap- 
ping potential y/ ad = r 2 /2. From Eq.Q, after inte- 
grating out the radial order parameter, the energy of the 
quasi- ID system is 



E 



/oo r 2 , 



N^uulipiflH 2 ] dz, (17) 



where Uu = 2au, u\ 2 — 2a\ 2 , and V i a (z) = a 2 z 2 /2 + 1. 
We analyze the ground state of the of repulsive TBEC 
aa,aij > trapped in quasi-lD traps in both miscible 
a i2 < \fd~\ai and immiscible a 12 > ^Ja\\a 22 domains. 
Without loss of generality we assume a 22 > an and 
hence, the first species, on account of lower repulsive 
mean field energy, has larger density at the center when 
Ni « N 2 . We consider 



,-( Z - 7 ) 2 /(26 2 ) 



(f 



(18) 



as our ansatz for the two order parameters with 
a, 7, b, /, c, 5, and d as variational parameters. Here 7 



4 



and S represent the location of the order parameter max- 
ima and center of mass motion. This is an apt ansatz as 
it describes the smooth transition between miscible and 
immiscible phases of the TBEC in a very natural way. 

In the miscible domain, the parameter c is a measure 
of flatness of the ip2 profile which arises from the larger 
intra-species repulsion energy. However, in the immisci- 
ble domain, it is the degree of separation between the two 
species due to the higher inter-species repulsion energy. 
The order parameters, like in previous subsection, satisfy 
the normalization conditions 



dz\^i{z)^ 



(19) 



where ipi(z) = \J a osc /Ni^>(za osc ), and hence the num- 
ber of independent variational parameters is reduced to 
five. From Eq.pjl), Eq.([l8|), and Eq.pjL we obtain the 



expression for energy as a function of five independent 
parameters. It is a complicated eighth degree polynomial 
(given in appendix) and can not be solved analytically. 
However, one can treat it as a nonlinear optimization 
problem and use numerical schemes like Nelder-Meade 
to find a solution [35]. Like in Eq.pl), the coupled ID 
GP equations 



3 = 1 



ipi = IMipi, (20) 



are obtained when £ = E — faN is variationally ex- 
trcmized with tp* as the variational parameters, where 
g u = N % Uii and = Nji 
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FIG. 2: Order parameters ipf of the two components from 
the variational calculations with equal intra-species scattering 
length of 51ao but no inter-species interaction (ay = 0). Solid 
(blue) and dashed (orange) correspond to the first and second 
species, respectively. The figure on the left hand side is for 
N\ — N2 = 1000, while that on the right hand side is for 
N\ — N2 — 100, 000. The trapping potential parameters are 
a — 0.02 and a osc = 9.566 x 10 _7 m. For (a) Niaii/a osc = 
2.821 and for (b) N t au/a osc = 282.137. 



1. Miscible domain with a\2 = 



two order parameters overlap and in the limiting case of 
a\2 = 0, the two maxima coincide if an = 022- Further- 
more, in the non-interacting limit ay — > 0, the profiles of 
the two order parameters are identical for the same num- 
ber of atoms (Ni = N2). For the special case of identical 
intra-species scattering lengths (an = 022), finite inter- 
species interaction (ai2 7^ 0), and identical number of 
atoms (N\ = N2), an appropriate ansatz is |36| 

M*) = Ae-^- B ^^ 2W ^e^ c+Dz+Ez2 \ 

ip 2 (z) = Ae-( z - B ^/^e l(c - Dz+Ez2 \ (21) 

originally introduced to describe dynamics of coupled 
solitons in nonlinear optical fibers [37]. Here, the pa- 
rameters A, B, C, D, E, and W are assumed to be time 
dependent and represent the amplitude, position, phase, 
wave number, chirp, and width of the Gaussian ansatz, 
respectively. In our ansatz too, a more symmetric choice 
of the order parameters 



Uz)^{h + c.^)e-^ 2 '^\ 



(22) 



can represent the special case mentioned here. The ex- 
pression of the energy is then much more complicated, 
and the situation considered being too restrictive, we 
do not consider this for further analysis and discussion. 
However, to examine the ground state geometry for a var- 
ied range of parameters, interaction strengths, and num- 
ber of atoms, the ansatz in Eq.(18) is an ideal choice. 

With our ansatz, when an = 022 ^ 0, ai2 = 0, and 
Ni = N 2 , the cz 2 dependence in ip 2 accounts for the 
self interaction. For this reason, the profile of ip2 is 
broader and more accurate, whereas tpi is a Gaussian 
and does not reflect the effect of the self interaction in 
the profile in an equally precise manner as i/V To ex- 
amine the density profiles as function of the nonlinearity, 
arising from the mean field interaction, consider TBEC 
with Ni = 100, 1000, and 10000. As a specific case take 
fl n = 0-22 = 51, which corresponds to 85 Rb [3J5], however, 
to begin with set 012 = 0. The later is experimentally not 
realizable in 85 Rb, but it is a good reference for a compar- 
ative study on the role of inter-species interactions. The 
order parameters for the different Ni with the chosen pa- 
rameters are as shown in Fig(2] As mentioned earlier, the 
i{>2 profiles are flatter and closer to the numerical values. 
Ideally, the two profiles should be identical as there is 
no inter-species interaction (ai2 = 0). So the difference 
between the profiles of the two order parameters is an 
indication of the error due to the Gaussian ansatz of ip\ . 
From the figures, the deviation grows as Ni is increased, 
which is expected as the mean field contribution to E is 
quadratic in N. In the weakly interacting domain c < 1 
and contributions to E from higher order terms of c are 
negligible. Retaining only the linear terms, the energy 
correction arising from the mean field is 



AE 



cfN 2 



25 2 - d 2 



da 2 (3d 4 + 12d 2 S 2 



As mentioned earlier, the inequality a 12 < \Ja\\a,22 de- 
fines the miscible domain. In this parameter range, the 



+4<5 4 ) + \/2df 2 5 22(d 2 + 4(5 2 ) 



(23) 



5 



For identical trapping potentials V\ — V 2 , the order pa- 
rameters ipi are centered at the origin ( z = 0) and hence 
5 ~ 0. The energy correction then simplifies to 



AE ' 



cdfN 2 



2j4 



(3o^d 



V2d 2 fg 22 - 1). (24) 



The three terms are the leading order corrections from 
the trapping potential, mean field, and kinetic energy, 
respectively. Considering that in quasi ID case a -C 1, 
the correction arising from c is a competition between 
the mean field and kinetic energy contributions; at lower 
values of 522, the kinetic energy correction dominates and 
AE is negative. As mentioned earlier for identical intra- 
species scattering lengths (an = 022), Eq.(pT|) is more 
appropriate ansatz and for parameters of Fig[2^b) the 
ground state geometry of the TBEC is shown in Figj3] 
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FIG. 3: Order parameters ipi of the two components from 
the variational calculations using Eq.(21l with equal intra- 



species scattering length of 51ao but no inter-species inter- 
action (aij = 0). Solid (turquoise) and dotted (black) cor- 
respond to the first and second species, respectively. The 
number of atoms and trapping potential parameters are same 
as those in Fig[2^b). 



the inter-species interaction energy is 

+b 8 (4f 2 +4cd 2 f + 3c 2 d 4 
+6cd 2 /) + 6 4 d 4 (24/ 2 + 3c 2 d 4 + 12cd 2 /) 



26 6 d 2 (8/ 2 + 3c 2 d 4 



• (25) 



Geometrically, with higher 012 the shape of i\) 2 near po- 
tential minima undergoes a smooth transition from Gaus- 
sian to flat Fig. |4][a-b) and then, to parabola with a lo- 
cal minima at the center of the trap Fig. [IJc-d). When 
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FIG. 4: Ground state wave functions for the two compo- 
nents of the TBEC with Ni = N 2 = 10000, 911 = 0.553, 
<?22 = 2gn, a = 0.02 and a osc — 9.566 x 10~ 7 m. Start- 
ing from left, the values of g\2 are 0.1(?22, 0.3(/22 in upper 
panel; and 0.5(?22 and 0.7<?22 in lower panel, respectively. The 
criterion, Njdij /a osc = gij/2 ~ 1, for the applicability of 
the variational ansatz is satisfied for each case. Turquoise 
(reddish-brown) and dark-green(orange) curves correspond to 
first and second component, respectively; solid curves are 
semi- analytic results, while dashed ones are the numerical so- 
lutions ofEq.pOl). 



2. Miscible domain with 0,12 > 

Finite inter-species interaction (a\ 2 7^ 0) modifies the 
density profiles of the two species in a dramatic way. As 
012 is increased, the species with the higher repulsion 
energy, second in the present case, is repelled from the 
center of trap. This is noticeable in the profiles of TBEC 
shown in Fig. Q . The lower density of the second species, 
around the center of the trap decreases inter-species den- 
sity product |V'i| 2 |' ( /'2| 2 & n d minimizes the total energy. 
However, this must be in proportion with the higher 
energy from trapping potential due to broader density 
profile. For identical trapping potentials (V~i = V 2 ), 
the profiles of the two species are centered at the origin 
(7 = 5 = 0) and symmetric. From Eq.(70) in appendix, 



-02 is flat around the center, to a very good approxima- 
tion ip' 2 (z) = ip 2 (z) = for z < 1. With our ansatz 
"02 (0) is not zero when a\ 2 > 0, so there are deviations 
from the actual solutions. Deviations are negligible when 
< 012 <ti 1, however, these grow and become quite 
prominent when tp 2 is constant around the center of the 
trap. 

At the center of the trap tp\ has the global maximum 
and V'i(O) = 0. Close to the center the order parameter 



Vh(0) + y<(0), 



(26) 



when z < 1. Gaussian ansatz of ipi defined in Eq.(18) 
then gives 



2 

z z a 



(27) 



() 



For the parameter domain where tp2 is constant close to 
the center for z < 1, we can write 



M*) « ^2(0). 



(28) 



After substituting these expressions of ipi in Eq.(20), the 



chemical potential of the two species obtained from the 
zeroth order equations are 



Mi = + 9na 2 +512^2(0) + 1, 



/'2 



522^2(0) +5210 2 



1. 



(29) 



Since ^(0) = 0, the first order terms are zero, and hence 
second order equation of ipi is 



3.9i 



-Vi(oM'(Q) 



^(o)+f*(oK(o) = ^i'(o). 

(30) 
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Further simplifications provide the condition 
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5ii = 



1 



2a 2 V 26 2 
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(31) 



on the intra-species interaction of the first species. Sim- 
ilarly, the second order equation of ip2 is 



a 

Y 



■<? 2 iVi(0)Vi(0) 



Mo) = 0. 



This defines the inter-species interaction as 



.921 



a 2 b 2 
2a 2 ' 



(32) 



(33) 



to obtain a constant profile of -02 around center of the 
trap. In other words when this condition holds true, the 
repulsive inter-species interaction balances the effect of 
confining potential. And the net outcome is an effective 
potential around the center which is constant. The ip2 
profile acquires a local minimum at the center when 321 > 
a 2 b 2 1 (2a 2 ), our ansatz is then the appropriate form of ip2- 



3. Immiscible domain: symmetric profile 



With harmonic trapping potentials the densities, to 
begin with, are not homogeneous in the miscible phase. 
Even when a 12 < *JawOyi, miscible domain discussed 
previously, there is a separation between the two maxima 
of the densities. When N\ ~ JV2 and an » 022, the 
density profiles are symmetric about the center (5 = 7 — 
0) and total energy of the TBEC is 
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(34) 



where E^ 1 is the inter-species interaction energy de- 



fined in Eq.(25). For the TBEC in square well de- 



scribed earlier, the criterion for phase separation |11) is 
the emergence of a dip in the total density pi{z)+p 2 (z) — 
\tpi(z)\ 2 + \i/j2{z)\ 2 at the interface. Such a criterion does 
not apply for TBEC in harmonic potentials. Instead we 
define phase separation as the state when the maxima of 
the densities are well separated. To quantify, define z e 
as the point where the densities of the two species are 
equal, which implies 



ae 



,/(26 ) 
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czl\e-<l^ d \ 



(35) 



For large a 12 , the order parameter ip2(z) is bimodal. The 
two maxima are solutions of V4( z ) = an d located at 



± 



^4cd 2 - / 



1/2 



(36) 



The maxima of %p2 are well separated from ipi when 
Zm > %e- We define this as the criterion for phase separa- 
tion of TBEC in harmonic potentials. This definition of 
phase separation is an appropriate one in very weakly in- 
teracting TBECs Na/a osc ~ 1, but not so in TF regime 
Na/a osc ^ 1, where / = is a better criterion to define 
phase separation. 



Phase separation occurs when the TBEC is immiscible. 
The density profiles of the two species are then spatially 
separated. Condition for immiscibility is 012 > ^Ja\\a22 
when the TBEC is confined within a square well poten- 
tial |llj . In which case, the trapping potential is flat 
and except for the surface effects, limited to within heal- 
ing length £ = ^/fi, 2 /(2m/i) from the boundary, there are 
no trapping potential induced density variations. The 
TBEC is then homogeneous in the miscible domain and 
phase separated in the immiscible domain, where the in- 
terface acquires a geometry which minimizes the surface 
energy. Here, the density distributions in the bulk is en- 
tirely dependent on the strength of interactions. And 
transition into phase separated domain is well defined. 



4- Imiscible domain: non- symmetric profile 

In imiscible domain, if one or both the intra-species 
non-linearities ga are sufficiently small, the stationary 
state geometry with symmetric ipi can have higher energy 
than the asymmetric stationary state with the two com- 
ponents lying side by side. There are two reasons for the 
emergence of asymmetric ground state geometries. First, 
this geometry has only one interface region, whereas the 
corresponding symmetric profile has two; as a result the 
interface energy is lower. And second, smaller kinetic en- 
ergy as ip[{z) is prominent at the interface and edges of 
the condensates. 
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FIG. 5: Ground state wave functions for the two compo- 
nents of the TBEC with TVi = iV 2 = 10000, gu = 0.553, 
<?22 = 2<7n, a — 0.02 and a osc = 9.566 x 10 _7 m. Starting 
from left, the values of gi2 are 0.9<?22, l-2<?22 in upper panel; 
and 1.4t/22 and 1.6322 in lower panel, respectively. The crite- 
rion, Njdij /dose = Sij/2 ~ 1, for the applicability of the vari- 
ational ansatz is satisfied for each case. Turquoise (reddish- 
brown) and dark-green (orange) curves correspond to first 
and second component, respectively; solid curves are semi- 
analytic results, while dashed ones are the numerical solutions 
of Eq.pOK 



With our ansatz, for the asymmetric profile c ~ 0, but 
S and 7 are nonzero. The total energy of the TBEC is 
then 



E 



a 2 ^ , f 2 N 2 



+ a 2 ba 2 N 1 {b 2 + 2 7 2 ) 



b d 

-df 2 a 2 N 2 {d 2 + 25 2 ) + V2a%N l9ll + V2df i N 2 g 22 
Aa 2 bde~^+ s f/^+^ 



(6 2 + d 2 ) 9 / 2 
+4b 6 d 2 f 2 +6b i d i f 2 )N l9l2 



(d 8 / 2 + Ab z d b f 2 + 6 s / 



(37) 



It must be emphasized that, with our present ansatz the 
asymmetric density profiles emerges very naturally as a 
function of a\ 2 . That is not the case with TF based semi- 
analytic methods which can not account for the existence 
of of asymmetric ground state geometries [5]. 

Typical ground state geometries of the TBECs ob- 
tained by this variational scheme are shown in Fig.([5| 
along with the respective non-linearities and trapping 
potential parameters mentioned in the caption of each 
figure. We observe that the present variational scheme 
also explains the existence of the asymmetric states, as 
is shown in Fig. ([5]), as the ground state geometry. 



C. Gram-Charlier expansion of \ipi\ 2 

As mentioned earlier, for the ideal case when all the 
non-linearities are small and equal, c is sj 0. In this 



limit, both the species have Gaussian distribution. In 
order to quantify the departure, brought about by inter- 
actions among the bosons, of the density distributions 
from the normalized Gaussian, we resort to an analysis 
based on Gram-Charlier series A. It is the expansion of a 
probability density function [38 in terms of the normal 
distribution. It can be used to analyze the departure 
from the normal distribution. If F(z) is a nearly nor- 
mal distribution with cumulants n r , then F(z) can be 
expressed as a series consisting of product of Hermite 
polynomials H^z) and K r . The truncated expression for 
Gram-Charlier A series up to fourth cumulant is 



F{z) 



-Z 2 /2 



24o- 5 V27r 
(z 4 - 6z 2 + 3) k 4 



24cr 4 + Aaz (z 2 - 3) k 3 + 



(38) 



where z — (z — n)/cr, with /1 and a as the mean and 
standard deviation of normal distribution. The departure 
of a distribution from normal are measured in terms of 
skewness ( 71 ) and kurtosis ( 72 )• These quantify the 
asymmetry and peakedness of the distribution function, 
and are defined as 



7i 



_K3_ 

3/2 • 



72 



K4 



(39) 



where second cumulant k 2 is equal to the variance a 2 . 
The [i and a 2 of the normal distribution are calculated 
from the numerically calculated ipi(z) using the defini- 
tions 

fjb = J z\^p(z)\ 2 dz, 

(z- Li) 2 \ip{z)\ 2 dz. 



For example, Fig. ^ shows the results of Gram-Charlier 
analysis of density profile for the second component of 
the TBEC in miscible and immiscible domains. As is 
evident from the figure, in the miscible domain, den- 
sity profile (solid-orange curve) is symmetric and rela- 
tively flatter (platykurtic) compared to normal distri- 
bution (dot-dashed reddish-brown curve) about origin, 
and consequently it has zero skewness and negative kur- 
tosis; whereas in second figure, the density distribution 
has a longer right tail (quantified by positive skewness) 
and narrow distribution about the mean (leptokurtic) as 
compared to normal distribution. From Gram-Charlier 
analysis, we can extract the coefficient of z A (excluding 
exponential part) and compare it with semi-analytic val- 
ues of c 2 . For both the cases in Fig.(|6|, the two values 
are same up to second decimal place. 
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FIG. 6: The solid orange curve is the approximate IV^I 2 
of the TBEC with Ni = 7V 2 = 10, 000, a osc = 9.566 x 



10" 



gn = 0.553, and 522 = 2gn, obtained by fitting 



Gram-Charlier A series over numerically obtained density 
profile (dashed dark-green curve). In (a)t/i2 = 0.6(722, while 
in (b) gi2 = 1.2g22- The dot-dashed reddish brown shows 
the density profile of normal distribution with respect to 
which 71 and 72 of numerical density profile have been cal- 
culated. For first figure, the variational parameters /, c, 
5, and d are 0.177, 0.001 ,0 , 8.675, respectively; while 
Gram-Charlier values are 0.185, 0, 0, 13.2836, respectively. 
For second figure, semi-analytic values of /, c, S, and d are 
0.267, 0, 8.529, 9.567,respectively; while Gram-Charlier val- 
ues are 0.270, 0, 7.358, 8.662, respectively. 



nth lattice site. Using the tight binding ansatz, energy 
functional of a TBEC in optical lattice is 

E = Xlfe [-2ifiViCiHCi(n + l) + e(n)iV i |C i (n)| 2 4 



■-NiUn)\ 4 



+ A 12 AT 1 |Ci(n)| 2 |C 2 (n)| 2 , (43) 



where A^, A12, K and e(n) are denned as 



A»t = 2a v iNi \ (f> n dz, 



Ai 2 = 2a 12 N 2 J (fadz, 

V z <^>„ • V z n+ i + <j) n V^{z)cj) n+1 



K 

e(n) 



1, 

2 

(V z 0„) 2 + V?(z) 



dz. 



dz, 
(44) 



To examine the stationary state of the system we consider 
Ci(n) is real while deriving above relation. 



III. TBECS IN OPTICAL LATTICES 

In this section, we analyze the ground state geometry 
of TBECs in optical lattices with the ansatz we have in- 
troduced. For this, we consider the optical lattice poten- 
tial, generated by a pair of orthogonally polarized counter 
propagating laser beams along axial direction, in pres- 
ence of the axisymmetric harmonic trapping potential. 
The period of the lattice potential is half of the laser 
wavelength. The net external potential experienced by 
TBEC (in scaled units) is the sum total of the two po- 
tentials, 



Vi(r,z) = ^(r 2 + a 2 z 2 ) + V Q cos 2 



2ttz 



(40) 



where Vq — sE r is the depth of potential well at each 
lattice site and A is the wavelength of the laser. Here, 
E r = (27r/A) 2 /2 is the recoil energy of laser light photon 
and s is the lattice depth scaling parameter. As in the 
previous section, in weakly interacting regime, the TBEC 
in cigar shaped traps is like a quasi-lD system and its 
energy is given by Eq.(17). The axial trapping potential 



V cos 



2ttz 



In tight binding approximation [33 



(41) 



(42) 



where n is index of lattice sites, 4> n { z ) — 4>{z— nX/2— A/4) 
is the wave function with amplitude C( n ) localized at 



A. Ground state 

In the weakly interacting regime, the localized wave 
function <fi n can be taken as ground state wave function 
of the lattice site, 



4>n{z) 



1 



= -( Z -nA/2-A/4)/(2 ( T 2 ) 



(45) 



where a = vA/ (8tt 2 Vq) 1 ^ 4 . Ideally, the energy associated 
with each individual lattice sites can be minimized to 
evaluate cj) n (z). Neglecting the energy due to harmonic 
trapping potential and interaction energy in comparison 
to potential energy due to optical lattice potential, the 
energy of the nth lattice site is approximately 




Vn 



2 r 



x-~(2n+l) 



dz, 



(46) 

where we have used Vq cos 2 {2itz/X) s» V (2tt/X) 2 (z — 
nX/2 — A/4) 2 as the potential at the nth lattice site. The 
expressions for e(n) and K are 



e(n) 



K 



Vn 



l-e 



1 



2 V" J ' 2,2a 2 

+a 2 cr 2 \ 2 (l + 2n) 2 ] , 

^ (V 



8a 2 a 4 



2 

2 \2 



1 + e 

x(l + n^-A 2 + 8(f7 2 + a 2 o- 6 



1 

32ct 4 



(47) 
[4a 4 A 2 a 2 

(48) 



Similarly, the energy per boson E = E/(Ni + N 2 ) can be 
calculated, however, as to be expected the expression is 
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FIG. 7: The stationary state solution of the TBEC with JVi = 
N 2 = 10,000, g n = 0.465, g 22 = 2gn, gv> = ffai = 0.8 fl2 2, 
a = 9/92, a oac = 1.137 x 10" 6 m, A = 0.7, and s = 6.0 as the 
non-linearity and trapping potential parameters, respectively. 
The criterion, Njdij /a oac = gij/2 ~ 1, for the applicability of 
the variational ansatz is clearly satisfied. The left and right 
figures are, respectively, the profiles obtained by variational 
minimization and numerically solving quasi- ID GP equation. 



FIG. 8: The stationary state solution of the TBEC with JVi = 
N 2 = 5,000, = 0.233, g 22 = 2gn, gi 2 = g 2 i = 2g 22 , 
a = 9/92, a osc = 1.137 x 10~ B m, A = 0.7, and s = 6.0 as the 
non-linearity and trapping potential parameters, respectively. 
The criterion, Njdij /a osc = gij/2 ~ 1, for the applicability of 
the variational ansatz is clearly satisfied. The left and right 
figures are, respectively, the profiles obtained by variational 
minimization and numerically solving quasi-lD GP equation. 



complicated and long. Following our ansatz, we consider 

&(n) = ae-(""^ 2 /(26 2 ) ) 

C 2 (n) = (f + cn 2 )e-^ 2 ^ d2 \ (49) 

as the envelope profiles with five independent parame- 
ters. The length scale associated with t he lattice poten- 
tial, in scaled units, is ai att = U uj /u)i , where u>i is the 
laser frequency of the optical lattice. Considering that 
the frequency of the harmonic trapping potential u> is at 
the most Khz and u>i is in the optical region autt <C 1. 
Hence, there are large number of lattice sites within the 
envelope profile. This implies that b 3> cr and d 3> er, 
where a is the width of the lattice ground state <fi n . We 
can replace summation over lattice sites with integration 
over the same variable in Eq.(43|. In fact, due to very 
weak trapping along axial direction, b and d are of the 

order of a few a osc even for very weak non-linearities. It 

2.2 2 .2 

implies that e" 71 " — > and e" 71 " — > 0, and in this 

limit the error incurred in replacing summation by in- 
tegration approaches zero. From Eq.(43), Eq.(44|, and 



Eq. ( 49 1 , we calculate the energy per particle which can be 



minimized to determine the variational parameters. The 
typical variational results for equal number of atoms of 



two species are shown in Fig. (11 1 



Symmetric and N on- symmetric Profiles 

Similar to TBECs in harmonic trapping potentials, 
ground sate geometry of binary condensates in optical 
lattices can be symmetric or asymmetric. In miscible 
and weakly segregated domain ( coherence length of the 
first component is smaller than the penetration depth of 
the second component [TT] ), the ground state geome- 
tries are symmetric. For example, a typical symmetric 
ground state geometry of TBEC in quasi-lD optical lat- 
tice is shown in Fig.Q. 

In strongly segregated domain ( coherence length of 
the first component is greater than the penetration depth 



of the second component), the ground state geometries 
can be asymmetric. A typical asymmetric ground state 
geometry of TBEC in quasi-lD optical lattice is shown 
in Fig.®. 



B. Profile narrowing with lattice depth 

The density profiles of the two components in TBEC is 
more compact with deeper lattice potential. To analyze 
this, we consider an ideal Bose gas ( non-interacting ) 
trapped in a triple ID potential well superimposed with 
a weak harmonic potential 



V(x) = < 



1 2 2 

-a x 
2 

1 2 2 

-ax 
2 

1 2 2 

-a x 
2 

00 



+ l -{5\x + df 



if 



+ 2^ 
+ 2/3 2 (x- 



if 



3d 

— x < 
2 

d d 
2 2' 



d 
2' 



d) 2 if - < x ^ 



3d 

y 



if \x\ > 3d/2, 



where a and are the parameters of harmonic and pe- 
riodic potentials, respectively, and d is the spatial extent 
of each well. This is most basic model to examine the 
nearest neighbour effects in a lattice. The central well, 
in particular, which has two nearest neighbours is a rep- 
resentation of lattice sites in optical lattices. In the ex- 
periment set up used by Cataliotti et. al. [30], the trap 
geometry is quasi- ID with Lu z /uj r — 9/92. For further 
analysis the equations are scaled in terms of the trans- 
verse oscillator length of the harmonic potential realized 
in aforementioned experimental set up, which is equiv- 
alent to setting a — 9/92. Accordingly, j3 and all the 
coordinates, from here after, are in scaled units. The 
potential considered is a simplified model and disconti- 
nuities are present at the boundary of two neighbouring 
wells. However, the ground state energy is much lower 
than the barrier height and it describes the underlying 
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physics very well. In the tight-binding approximation, 
the total wave function is 



(51) 



where i = —1,0, and 1 represent the left, central, and 
right well, respectively. Here the normalization is 



\ip(x)\ 2 dx = 1. 



(52) 



We consider f3 > a, which follows from the tight bind- 
ing approximation. The width of each of these localized 
wave functions are calculated by minimizing the localized 
energies of each well 



E; 



Ni 



|V&| 2 + V5(aO|& 



dx, 



where, from the previous definition 

Vi(x) = ^x 2 + ^/3 2 (id + x) 2 . 



(53) 



(54) 



The range of each well, as defined earlier, is (2i — l)d/2 ^ 
x ^ (2i + l)d/2. Neglecting the nearest neighbour over- 
laps of wave the functions, the integration limits can be 
considered as — oo to +oo. The energies are then 



Eo = 2 V, 
E±i = Eq + SE 



(55) 



where r\ = \J a 2 + f3 2 and 



2 J2 



SE 



a 2 d 



(56) 



is the energy difference between side and central wells 

2.5 [ 




200 220 240 260 280 300 320 

FIG. 9: Variation of v as a function of j3. 

For the combined system of the three wells, in tight 
binding approximation, we only consider tunneling be- 
tween adjacent wells. The lowest eigen energy is the 



chemical potential of the system and the eigen vector cor- 
respond to the probability amplitudes for the occupancy 
of each well. As defined earlier, E_i, E , and Ei are the 
energies localized in left, central, and right well, respec- 
tively. Define the tunneling matrix element between left 
and central well as 



V0_i ■ V</»o + <f>-i V{x)cj> 



dx, 



and similarly, for the central and right wells as 
T. 



-V0 O • V0i + 0o^)0i 



dx. 



(57) 



(58) 



For the symmetric case considered here v — w, and after 
evaluating the integrals 



16y/nrj 



-id 



1 1 - 2e 



3d 2 ?; 



- 2e 



4d~ t? 



- 1 + 2e 



3d 2 17 



'nrje 



■nrje 



4d 2 ri 



[27 ? + d 2 (5/3 2 -3)] Erf(d^) 
[2?y - rf 2 (ll/3 2 + 3)] Erf (2dy/ff) 



(59) 



In the above expression Erf(---) represents error func- 
tion. The Hamiltonian matrix of the system is then 



'E-i v 
H= \ v Eo v 
v Ei 



(60) 



The ground state of the system is the lowest energy eigen 
vector obtained from diagonalizing the Hamiltonian ma- 
trix. This is equivalent to solving the secular equation 
and amounts to calculating the roots of a cubic polyno- 
mial, which is possible analytically. The eigen values, 
in increasing order of magnitude, of the Hamiltonian in 
Eq.((60l) are 



-(E0+E1 — E12), Ei, -(E + Ex + E12); 



(61) 



where Ev> + 8v 2 . The eigen vector of the lowest 

eigen value is 



a±i 



a = 



V2v 


^8v 2 + 


SE {SE + 


E12) 


ys« 2 


+ SE (SE 


+ E 12 ) 



E z + 8v 2 



and the ratio of the probability amplitude of occupancy 
of central to side well is 



|o | _ SE(SE + E 12 ) + 8v 2 
\ a ±i\ 2v\/SE 2 + 8v 2 



(62) 



This energy difference SE is quite small for small values 
of d. For example based on Ref.[3(Tl, consider d = 0.35, 
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and P = 197.824. The associated energy difference and 
tunneling matrix element are 5E = 0.00049494 and v — 
2.18549, respectively. Hence, for SE <§C v, using Taylor 
series expansion Eq. ( 62 1 simplifies to 



•lAx) 



Igol 
l a ±i| 



VI- 



SE 
2u ' 



(63) 



Since 5E is independent of f3 and v decreases as f3 is in- 
creased (see Fig. ([9])), the above ratio is larger with higher 
values of f3. Using Eq.(63), Fig. (10) shows the variation 



of ratio of probability of occupancy of central to side 
wells with respect to v. This relation is not valid for 
very large values of /? when SE v does not hold true. 
For trapping potentials without the harmonic component 
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FIG. 11: The plots show the ground state wave function 4>(x) 
(solid curve) and probability den sity ^(x)! 2 (dashed curve) 
for potential considered in Eq.( 50 1 For figure on left side, a = 



9/92, d = 0.35 and /3 = 197.82 (s = 6), while for figure on 
right side, a = 9/92, d = 0.35 and P = 312.79 (s = 15) 
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have occupancy probability of 0.249976, slightly lower 
probability of occupancy than a = case. The trend 
continues with further increase of (3, and for f3 = 1000 
there are almost no atoms in the potential wells at the 
wings. Thus with increase in value of /3, atoms migrate 
from the wings to central well and lead to the squeezing of 
the condensate. The squeezing of the condensate profile 



is evident from the comparative study of plots in Fig( 11 1 



In three level approximation, a triple well problem can 
also be analyzed by considering the following classical 
Hamiltonian l4"2l 



FIG. 10: The plot shows the po/p±i as a function of v for 
the values of (3 ranging from 197.82 to 312.79. In terms of the 
ratio of barrier height to recoil energy of laser light photon 
s used in Ref . |30| . the above mentioned range of /3 values is 
equivalent to s ranging from 6 to 15. 

(a = 0) , purely lattice potential, the three eigenvalues are 
identical when the hopping across adjacent wells is zero 
(v = w = 0). This solution correspond to all the atoms 
confined to one of the potential wells. At this stage for 
further reference we define the occupancy probability of 
the ith well as 



Pi 



<f)*ipdx 



(64) 



When v and w are non-zero, then the ground state has 
nonzero occupancies for all the three wells. The central 
well has probability of occupancy equal to half (p = 0.5) 
while the wells at the wings have probability equal to 
one fourth each (p±i = 0.25). More importantly, these 
values are are independent of (3. That is, no squeezing of 
the density profile occurs without the harmonic potential 
(a = 0). ' 

In the presence of a weak harmonic potential, there is 
an increase in p when /3 is increased. For = 197.82, 
equivalent to s = 6 in Ref. [30] , the occupancy probability 
of the central well is 0.500047, this is marginally larger 
than the a = case. On the other hand, the wings 



n 



1 



g [n 2 _ 1 + (1 - n_i 



ni, 



e_i + niei 



— n-i — ni (vy/n-i cos 6>_i + w^Jn\ cos 9\) ,(65) 
; \aj\ 2 and 6j = arg(a2) — arg(cij) are canonical 



where rij = 

conjugate variables satisfying following equations 



drij 
lit 



and 



09 



dt 



on 

Orij 



(66) 



and e_i and e± are the zero point energies of potential 
wells at the flanks measured with respect to zero point 
energy of central well. The stationary states of the sys- 
tem can be calculated by by solving 



drij 
~dT 



dt 



(67) 



In the absence of the harmonic potential e_i and ei 
are zero and change in tunneling amplitude has no ef- 
fect on the occupancy of three wells. But the situa- 
tion is dramatically different with the harmonic poten- 
tial, e_i = t\ 7^ 0. In this case probability of occupancy 
of central well increases with decrease in v. This is ex- 
perimentally realizable by increasing depth of lattice po- 
tential. Similarly for four potential wells, the occupancy 
of the two central wells (whose occupancies are equal) 
grows as we increase /3 for non zero a. 
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IV. CONCLUSIONS 

We have studied the stationary state properties of 
the TBECs both in miscible and imiscible domains us- 
ing a variational ansatz. The stationary state geome- 
tries obtained using the variational ansatz are in very 
good agreement with numerical results, especially for 
very weakly interacting TBECs, i.e., Na/a osc is of the 
order 1. We have also quantified the departure of the 
variational ansatz based semi-analytic results from the 
numerical ones using Gram- Char lier analysis. Besides 
harmonic trapping potentials, the present ansatz can also 
be used to study the TBECs in deep optical lattices where 
coupled CP equations can be mapped into coupled dis- 
crete non-linear Schrddinger equations. In optical lat- 



tices, the density profiles of the component species of 
the TBEC get squeezed with the increase in the height 
of the potential barrier between adjacent wells. We have 
explained this phenomenon using a very simple triple well 
potential model. 
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APPENDIX 



r„2 



+ a 2 ba 2 (b 2 + 7 2 ) + V2a 4 b gil 



(68) 



^{^f-^-2S 2 ) + c 2 (7d 4 



+20d 2 5 2 + A5 4 

rl2 x2 



E 



12 



64d[4/ 2 (rf 2 + 25 2 ) +4c/(3rf 4 
+12d 2 S 2 + AS 4 ) + c 2 (15d 6 + 90d 4 5 2 + 60d 2 S 4 
+85 e )]a 2 + V2d[256.f 4 + 256c f 3 (d 2 + AS 2 ) 
+96c 2 / 2 (3d 4 + 24d 2 <5 2 + 165 4 ) + I6c 3 f(l5d e 

+180d 4 S 2 + 2A0d 2 S 4 + 64<5 6 ) + c 4 (105d 8 + 1680d 6 (5 2 



/Tra 2 bdNigi 2 e » 2 + d2 



{Ad\f + C1 2 f 



A(b 2 + d 2 ) 9 / 2 
+Ab 2 d 6 [4/ 2 + cf(d 2 + 4 7 2 - 4 7 <5) + c 2 7 2 (3d 2 
-4ryS)] + b 8 [Af 2 + Acf(d 2 + 25 2 ) + c 2 {U 4 + I2d 2 6 2 
+A5 4 )] + 26 6 d 2 [8/ 2 + c 2 (3d 4 - 8 7 <S 3 - \2d 2 b 1 
+6d 2 S 2 ) + 2cf(3d 2 - 4(5 7 + AS 2 )] + b 4 d 4 [24/ 2 
+3c 2 (d 4 + 4d 2 7 2 - 8d 2 ~/5 + 8 1 2 5 2 ) + 4c/(3d 2 
+2 7 2 - 8 7 <5 + 25 2 )} } (70) 



3d 5 + 12d 3 6 2 + AdS 4 



-d 3 f - 2dfS 2 



(3d 5 y^ 



2d 6 f 2 n + 12d 3 y^S 2 - 8d 4 fir5 2 + Ad^S 4 ) 



is 1/2 



(71) 
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